TopographicIndex Subroutine

public subroutine TopographicIndex(dem, fdir, facc, xcoord, ycoord, tivalue, tigrid)

Topographic index TI takes into account both a local slope geometry and site location in the landscape combining data on slope steepness S and specific catchment area As:

   TI = ln (As/S)

For grid structures the specific catchment area can be obtained by dividing the contributing area of the cell by the effective contour length. This is the length of contour line within the grid cell over which flow can pass. It is calculated as the length of the line through the grid cell centre and perpendicular to the aspect direction.

References:

P. J. J. DESMET & G. GOVERS (1996): Comparison of routing algorithms for digital elevation models and their implications for predicting ephemeral gullies, International Journal of Geographical Information Systems, 10:3, 311-331

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(in) :: dem

digital elevation model
type(grid_integer), intent(in) :: fdir

flow direction
type(grid_real), intent(in), optional :: facc

flow accumulation in cell unit
real(kind=float), intent(in), optional :: xcoord

x coordinate of cell for computing index
real(kind=float), intent(in), optional :: ycoord

y coordinate of cell for computing index
real(kind=float), intent(out), optional :: tivalue

index of cell with coordinates (x,y)
type(grid_real), intent(out), optional :: tigrid

topographic index

Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: alpha

local aspect (rad)
type(grid_real), public :: area

contributing area (m^2)
real(kind=float), public :: as

specific catchment area (m^2 / m)
type(grid_real), public :: aspect

aspect (rad)
real(kind=float), public :: cellsize

cell size (m)
integer(kind=short), public :: i
integer(kind=short), public :: j
type(grid_real), public :: slope

slope (rad)
real(kind=float), public :: x

factor converting cellsize along direction... ...perpendicular to the aspect direction

Source Code

SUBROUTINE TopographicIndex &
!
(dem, fdir, facc, xcoord, ycoord, tivalue, tigrid)

IMPLICIT NONE

!Argumentd with intent(in):
TYPE (grid_real), INTENT(IN) :: dem !!digital elevation model
TYPE (grid_integer), INTENT(IN) :: fdir !!flow direction
TYPE (grid_real), OPTIONAL, INTENT(IN) :: facc !!flow accumulation in cell unit
REAL (KIND = float) , OPTIONAL, INTENT(IN) :: xcoord !! x coordinate of cell for computing index
REAL (KIND = float) , OPTIONAL, INTENT(IN) :: ycoord !! y coordinate of cell for computing index

!Arguments with intent(out):
REAL (KIND = float) , OPTIONAL, INTENT(OUT) :: tivalue !! index of cell with coordinates (x,y)
TYPE (grid_real), OPTIONAL, INTENT(OUT) :: tigrid !!topographic index

!local declarations
TYPE (grid_real) :: area !! contributing area (m^2)
TYPE (grid_real) :: aspect !!aspect (rad)
TYPE (grid_real) :: slope !!slope (rad)
REAL (KIND = float) :: as !! specific catchment area (m^2 / m)
REAL (KIND = float) :: x !!factor converting cellsize along direction...
                         !!...perpendicular to the aspect direction
REAL (KIND = float) :: alpha !! local aspect (rad)
REAL (KIND = float) :: cellsize !!cell size (m)
INTEGER (KIND = short) :: i, j
!------------------------------end of declaration -----------------------------

!allocate topographic index grid
CALL NewGrid (tigrid, dem)

IF ( .NOT. PRESENT(facc) ) THEN
   !compute contributing area
   CALL FlowAccumulation (fdir, area)
ELSE
    CALL  NewGrid (area, fdir)
    DO i = 1, area % idim
      DO j = 1, area % jdim
          IF (area % mat (i,j) /= area % nodata ) THEN
              area % mat (i,j) = facc % mat (i,j) ! (* CellArea (facc, i, j))
          END IF
      END DO
    END DO
END IF

!compute aspect
CALL DeriveAspect (dem, aspect)

!compute slope
CALL DeriveSlope (dem, slope)

!if present tivalue compute topographic index only for cell with coordinate xcoord, ycoord
IF (PRESENT (tivalue) ) THEN
    !check wether coordinates of cell is given
    IF ( PRESENT (xcoord) .AND. PRESENT (ycoord) ) THEN
        CALL GetIJ (xcoord, ycoord, fdir, i, j)
        
        alpha = aspect % mat (i,j)
        x = 1. / MAX(ABS(SIN(alpha)),ABS(COS(alpha)))
        cellsize = ( CellArea (fdir, i, j) ) ** 0.5 ![m]
        as = area % mat (i,j) / (cellsize * x)
        IF ( TAN (slope % mat (i,j)) == 0.) THEN
            tivalue = LOG (as / 0.0001 )
        ELSE
            tivalue = LOG (as / TAN (slope % mat (i,j)))
        END IF
 
    ELSE
       CALL Catch ('error', 'Morphology', 'missing coordinate for compuing topographic index of given cell')
    END IF
    
END IF

    
!If present tigrid, compute topographic index for every cell of the grid
IF (PRESENT (tigrid) ) THEN

        DO i = 1, tigrid % idim
          DO j = 1, tigrid %  jdim
            IF (tigrid % mat (i,j) /= tigrid % nodata) THEN
              alpha = aspect % mat (i,j)
              x = 1. / MAX(ABS(SIN(alpha)),ABS(COS(alpha)))
              cellsize = ( CellArea (tigrid, i, j) ) ** 0.5 ![m]
              as = area % mat (i,j) / (cellsize * x)
              IF ( TAN (slope % mat (i,j)) == 0.) THEN
                 tigrid % mat (i,j) = LOG (as / 0.0001)
              ELSE
                 tigrid % mat (i,j) = LOG (as / TAN (slope % mat (i,j)))
              END IF
     
            END IF 
          END DO
        END DO
        
END IF


!purge
CALL GridDestroy (area)
CALL GridDestroy (aspect)
CALL GridDestroy (slope)

RETURN

END SUBROUTINE TopographicIndex